Abstract
Hemos analizado con las herramientas proporcionadas en el curso de Modelos Estadísticos el conjunto de datos, Cheddar, distribuido en la librería Faraway de R. Para ello hemos utilizado diversas técnicas de regresión lineal y no lineal.En un estudio de queso Cheddar realizado en el Valle de Latrobe (Victoria, Australia), se estudiaron muestras de queso en las que se analizó su composición química y fueron dadas a probar a distintos sujetos para que valoraran su sabor. Los valores asignados a cada queso son el resultado de combinar las distintas valoraciones.
El DataFrame cheddar de la librería faraway consiste de 30 muestras de queso Cheddar en las que se ha medido el sabor (taste) y las concentraciones de ácido acético (Acetic), ácido sulfhídrico (H2S) y lactosa (Lactic).
Tenenemos un conjunto de datos en el que se recogen observaciones de una cata de quesos, nuestras variables son:
· Taste: una valoración subjetiva de los jueces.
· Acetic: la concentración de ácido acético en un queso de terminado en esca la logarítmica
· H2S: la concentración de sulfito de hidrógeno en escala logarítmica.
· Lactic: Concentración de ácido láctico
A lo largo del documento hacemos uso de las siguientes librerias de R:
library(faraway)
library(leaps)
library(MASS)
library(PASWR)
library(car)
library(ggplot2)
library(plyr)
library(GGally)
library(corrplot)
library(plotly)
library(scatterplot3d)
library(cowplot)
Vamos a utilizar el dataset Cheddar. Cargamos los datos y enseñamos las primeras observaciones.
# load cheddar cheddar
data(cheddar)
head(cheddar)
A continuación, comprobamos que no hay entradas vacias
# check if any entry has a missing value (NA)
any(is.na(cheddar))
## [1] FALSE
sapply(cheddar, class)
## taste Acetic H2S Lactic
## "numeric" "numeric" "numeric" "numeric"
Al ser todas numéricas (cuantitativas) no hay que hacer encoding a variables binarias. Con estas variables famos a intentar explicar cómo los valores observados de una variable Y (taste) dependen de los valores de otras variables (Acetic, H2S, Latcic), a través de una relación funcional del tipo Y = f(X) También vamos a intentar predecir el valor de la variable Y para valores no observados de las variables X.
Usamos el número de observaciones para determinar los conjuntos de train y test.
numObs <- dim(cheddar)[1]
numObs
## [1] 30
Tenemos 30 observaciones, ¿Podemos suponer que la distribución de las variables es normal?,¿Tenemos alguna en la que falten datos?, ¿Tenemos outliers?, etc.
Las anteriores preguntas las trataremos en las siguientes secciones.
# Para asegurar que sea reproducible
set.seed(1)
train <- sample(c(TRUE, FALSE),
size = nrow(cheddar),
replace = TRUE,
prob = c(0.7, 0.3))
test <- (!train)
Hacemos un pequeño estudio preliminar de nuestras variables:
plot(cheddar)
summary(cheddar) # Distribución de las variables
## taste Acetic H2S Lactic
## Min. : 0.70 Min. :4.477 Min. : 2.996 Min. :0.860
## 1st Qu.:13.55 1st Qu.:5.237 1st Qu.: 3.978 1st Qu.:1.250
## Median :20.95 Median :5.425 Median : 5.329 Median :1.450
## Mean :24.53 Mean :5.498 Mean : 5.942 Mean :1.442
## 3rd Qu.:36.70 3rd Qu.:5.883 3rd Qu.: 7.575 3rd Qu.:1.667
## Max. :57.20 Max. :6.458 Max. :10.199 Max. :2.010
Ploteamos las gráficas de dispersion entre la variable respuesta taste y las variables predictoras.
layout(matrix(1:3, nrow = 1))
# plot taste ~ Acetic
plot(Acetic, taste,
main = "Relación entre Taste y Acetic",
xlab = "Acetic", ylab = "Taste",
pch = 20, frame = FALSE)
# plot taste ~ H2S
plot(H2S, taste,
main = "Relación entre Taste y H2S",
xlab = "H2S", ylab = "Taste",
pch = 20, frame = FALSE)
# plot taste ~ Lactic
plot(Lactic, taste,
main = "Relaciónn entre Taste y Lactic",
xlab = "Lactic", ylab = "Taste",
pch = 19, frame = FALSE)
Podemos observar que la que aperentemente guarda una menor relación lineal
con taste es la variable Acetic, esto será comprobado con distintos tests.
Simplificamos nuestra notación para las variables, la que intentaremos predecir es taste.
Definimos nuestro modelo completo:
x <- model.matrix( ~ Acetic + H2S + Lactic, data = cheddar[train,])
betahat <- solve(crossprod(x, x), crossprod(x, taste))
betahat <- c(betahat)
betahat
## [1] -14.016212 -5.066965 3.217036 31.938449
model.all.lm <- lm(taste ~ ., data = cheddar[train,])
model.all.lm$coefficients
## (Intercept) Acetic H2S Lactic
## -14.016212 -5.066965 3.217036 31.938449
Que evidentemente son los mismos.
Estudiemos preliminarmente si es un modelo lineal adecuado, para ello comprobaremos las hipótesis estándar del modelo lineal de regresión
shapiro.test(resid(model.all.lm))
##
## Shapiro-Wilk normality test
##
## data: resid(model.all.lm)
## W = 0.97755, p-value = 0.8865
Observamos que estamos en la hipótesis de que el error nuestro modelo se distribuye de manera normal. (p-value = 0.8865 > 0.05)
¿Tienen todas las variables un impacto relevante en el modelo?
summary(model.all.lm)
##
## Call:
## lm(formula = taste ~ ., data = cheddar[train, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.7641 -4.0925 -0.7499 4.3790 17.7545
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -14.016 20.490 -0.684 0.5032
## Acetic -5.067 5.120 -0.990 0.3363
## H2S 3.217 1.472 2.185 0.0432 *
## Lactic 31.938 12.852 2.485 0.0237 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.109 on 17 degrees of freedom
## Multiple R-squared: 0.7251, Adjusted R-squared: 0.6765
## F-statistic: 14.94 on 3 and 17 DF, p-value: 5.103e-05
Obervamos el p-valor de A nos indica que con casi toda seguridad A no tiene impacto real en el modelo ( > 0.94)
Los p-valores son lo suficientemente bajos como para rechazar varianza constante
Una vez hemos concluido que aunque estamos en las hipótesis de regresión lineal el modelo completo a pesar de ser el más complejo probablemente da resultados similares a otro más simple.
corplot <- ggpairs(cheddar[train,], progress=FALSE)
corplot
mat_cor <- cor(cheddar, method = "pearson")
corrplot(mat_cor, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)
anova(model.all.lm)
Todos nuestros p-valores son adecuados a un nivel = 0.05
Para comprobarlo basta realizar el siguiente test:
outlierTest(model.all.lm)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 12 2.458893 0.025711 0.53993
Concluimos con un nivel \(\alpha\) = 0.05 que no tenemos ningún outlier en nuestra muestra.
Esto se puede comprobar graficamente a través del siguiente gráfico
influenceIndexPlot(model.all.lm)
Observamos que hay un pequeño grupo de observaciones que aunque no sean outliers pueden alterar la claidad del modelo.
Partimos del modelo completo estudiado en la sección anterior y aplicamos con = 0.05.
drop1(model.all.lm, test = "F")
Quitamos Acetic del modelo debido que su p-valor es > 0.05
model.update1 <- update(model.all.lm, . ~ . - Acetic)
drop1(model.update1, test = "F")
model.update2 <- update(model.all.lm, . ~ . - H2S)
drop1(model.update2, test = "F")
Todos los p-valores son menores que = 0.05 por lo tanto hemos concluido
model.backward.lm <- lm(taste~Lactic, data = cheddar[train,])
Modelo resultante: taste ~ Lactic.
SCOPE<-(~.+Acetic + H2S + Lactic)
model.inicial <- lm(taste~1,data= cheddar[train,])
add1(model.inicial,scope=SCOPE,test="F")
Actualizamos añadiendo Lactic
model.updatei1 <- update(model.inicial, .~. + Lactic)
add1(model.updatei1,scope=SCOPE,test="F")
Con nivel de significación = 0.05 este sería nuestro modelo final.
model.forward.lm<-lm(taste~Lactic, data= cheddar[train,])
Modelo resultante: taste ~ Lactic.
A este modelo, taste ~ Lactic le vamos a dar un nombre, posteriormente estudiaremos los problemas que presenta.
model.final1 <- lm(taste ~ Lactic, data = cheddar[train,])
En esta subsección trataremos de encontrar un candidato a mejor modelo, construyendo nuestros modelos usando distintos enfoques.
models <- regsubsets(taste ~ ., data = cheddar[train,])
summary(models)
## Subset selection object
## Call: regsubsets.formula(taste ~ ., data = cheddar[train, ])
## 3 Variables (and intercept)
## Forced in Forced out
## Acetic FALSE FALSE
## H2S FALSE FALSE
## Lactic FALSE FALSE
## 1 subsets of each size up to 3
## Selection Algorithm: exhaustive
## Acetic H2S Lactic
## 1 ( 1 ) " " " " "*"
## 2 ( 1 ) " " "*" "*"
## 3 ( 1 ) "*" "*" "*"
MR2adj <- summary(models)$adjr2
MR2adj
## [1] 0.6197312 0.6769081 0.6765346
which.max(MR2adj)
## [1] 2
summary(models)$which[which.max(MR2adj), ]
## (Intercept) Acetic H2S Lactic
## TRUE FALSE TRUE TRUE
Modelo resultante: taste ~ H2S + Lactic.
MCp <- summary(models)$cp
MCp
## [1] 5.336570 2.979216 4.000000
which.min(MCp)
## [1] 2
summary(models)$which[which.min(MCp), ]
## (Intercept) Acetic H2S Lactic
## TRUE FALSE TRUE TRUE
Modelo resultante: taste ~ H2S + Lactic.
MBIC <- summary(models)$bic
MBIC
## [1] -15.29253 -16.80519 -14.93673
which.min(MBIC)
## [1] 2
summary(models)$which[which.min(MBIC), ]
## (Intercept) Acetic H2S Lactic
## TRUE FALSE TRUE TRUE
Modelo resultante: taste ~ H2S + Lactic.
stepAIC(model.all.lm, scope = SCOPE, k = 2)
## Start: AIC=96.35
## taste ~ Acetic + H2S + Lactic
##
## Df Sum of Sq RSS AIC
## - Acetic 1 81.26 1492.0 95.530
## <none> 1410.7 96.354
## - H2S 1 396.17 1806.9 99.551
## - Lactic 1 512.50 1923.2 100.862
##
## Step: AIC=95.53
## taste ~ H2S + Lactic
##
## Df Sum of Sq RSS AIC
## <none> 1492.0 95.530
## + Acetic 1 81.26 1410.7 96.354
## - H2S 1 361.58 1853.5 98.087
## - Lactic 1 443.44 1935.4 98.994
##
## Call:
## lm(formula = taste ~ H2S + Lactic, data = cheddar[train, ])
##
## Coefficients:
## (Intercept) H2S Lactic
## -31.122 3.054 25.210
Modelo resultante: taste ~ H2S + Lactic.
Nótese que los modelos obtenidos por metodos de criterios coinciden.
model.crit <- lm(taste ~ H2S + Lactic, data=cheddar[train,])
anova(model.crit, model.all.lm)
Con un valor de 0.3363 no podemos decir que los modelos sean signficativamente diferentes con ningún nivel habitual, interpretamos que el segundo modelo aunque más simple, no es significativamente peor. #### Observación
A este , taste ~ H2S + Lactic le vamos a dar un nombre, en la siguiente sección lo estudiaremos con detalle.
model.final2 <- lm(taste ~ H2S + Lactic, data = cheddar[train,])
En esta sección estudiaremos si nuestros modelos cumplen las condiciones necesarias de un modelo de regresión lineal.
Nuestro enfoque consistirá en un analisis gráfico, acompañado de tests estadísticos en los casos en los que se aprecie una discrepancia notable.
En la sección 3 se toma un enfoque naïve a la hora de construir los modelos, ya que no hemos estudiado si hay observaciones influyentes, podríamos tener una muestra que no es la adecuada para el estudio de nuestros datos.
Un modelo de regresión lineal debe satisfacer las siguientes hipótesis con nivel de siginificación adecuado:
plot(model.final1)
residuos <- resid(model.final1)
Observamos que este modelo tiene algunas características poco lineales, ya que en el primer plot hay una aparente ausencia de relación linea, lo comprobamos estadisiticamente.
plot(model.final2)
residuos2 <- resid(model.final2)
Observamos que este modelo tiene algunas características poco lineales, ya que en el segundo plot nuestros datos llegan a dispersarse considerablemente,lo comprobamos estadisiticamente.
A fin de exponer con una mayor claridad los resultados y ahorrar espacio, recogemos los resultados de los tests en una tabla en la que observamos si se verifican las hipotesis del modelo de regrsión lineal.
library(knitr)
## Warning: package 'knitr' was built under R version 4.1.3
res.aov <- aov(model.final1)
res.aov2 <- aov(model.final2)
shap1 <- shapiro.test(residuos)$p.value
shap2 <- shapiro.test(residuos2)$p.value
t1 <- t.test(residuos, mu = 0, alternative = "two.sided")$p.value
t2 <- t.test(residuos2, mu = 0, alternative = "two.sided")$p.value
aov1L <- summary(res.aov)[[1]][["Pr(>F)"]][1]
aov2L <- summary(res.aov2)[[1]][["Pr(>F)"]][1]
aov2H <- summary(res.aov2)[[1]][["Pr(>F)"]][2]
dw1 <- durbinWatsonTest(model.final1)[[3]]#CORREGIR PROBLEMILLA; NO SON ESTOS LOS P
dw2 <- durbinWatsonTest(model.final2)[[3]]
df_hip <- data.frame("Distribución_normal" = c(shap1,shap2, 0.05,"Shapiro-Wilk"),
"Media_0" = c(t1,t2, 0.05,"t-Test"),
"Varianza_no_constante_H2S" = c("No participa en modelo",aov2H,0.05,"ANOVA"),
"Varianza_no_constante_L" = c(aov1L,aov2L,0.05,"ANOVA"),
"No_Autocorrelación" = c(dw1,dw2,0.05, "Durvin-Watson"))
rownames(df_hip) <- c("model.final1","model.final2","p-valor","Test usado")
df_hip
A pesar de nuestras suposiciones iniciales, los modelos satisfacen todas las hipótesis de un modelo de regresión lineal.
Es interesante observar que los residuos del segundo modelo se ajustán mejor a una normal y la variación de la autocorrelación al añadir un segundo predictor.
Anteriormente hemos tratado que nuestra muestra no contiene outliers. ¿Pasa lo mismo con las observaciones influyentes?
Realizamos un estudio de estos últimos a través de distintos criterios: #####4.2.1 Estuudio de model.final1
X <- model.matrix(~ H2S + Lactic, data = cheddar)
H <- X %*% solve(t(X) %*% X) %*% t(X)
hii <- diag(H)
hCV <- 2 * 3 / 30
sum(hii > hCV)
## [1] 1
which(hii > hCV) # 6
## 6
## 6
dffitsCV <- 2 * sqrt(3 / 30)
dffitsmodel <- dffits(model.final1)
sum(dffitsmodel > dffitsCV)
## [1] 2
dfbetaCV <- 2 / sqrt(30)
dfbetamodel <- dfbeta(model.final1)
dfbetamodel
## (Intercept) Lactic
## 1 10.98453709 -6.857931612
## 2 0.04813185 -0.204187174
## 3 -0.50663219 0.739444044
## 5 1.52600073 -0.930155631
## 8 3.53325599 -2.954743648
## 9 0.50645416 -0.258946998
## 10 0.37463452 -0.499299660
## 11 -0.36774193 0.342168554
## 12 -5.72680931 4.515468710
## 13 -2.19024714 1.309997368
## 14 1.89785939 -0.955345774
## 16 -1.19914765 1.035728744
## 19 1.18937074 -1.238273837
## 22 -0.34585424 0.164146758
## 23 0.67946353 -0.123569461
## 24 -5.26050279 4.025271578
## 25 0.06574074 -0.032522575
## 26 -0.57342391 0.006350913
## 27 1.36540211 -1.203660261
## 28 -3.13839299 1.685507961
## 30 -1.99398262 1.070889972
sum(dfbetamodel[, 1] > dfbetaCV)
## [1] 9
sum(dfbetamodel[, 2] > dfbetaCV)
## [1] 7
#sum(dfbetamodel[, 3] > dfbetaCV) # Yo diria que esta se quita para evitar error
which(dfbetamodel[, 1] > dfbetaCV)
## 1 5 8 9 10 14 19 23 27
## 1 4 5 6 7 11 13 15 19
#which(dfbetamodel[, 3] > dfbetaCV) # Yo diria que esta se quita para evitar error
Una vez hemos localizado estas observaciones, en este caso las hay, las retiramos de la muestra con el fin de tener una muestra adecuada para trabjar
obs.out <- c(6, 7, 8, 12, 15)
cheese <- cheddar[-obs.out, ]
Ahora ya estamos en condiciones de realizar un estudio adecuado de las condiciones de regresión lineal.
model.exh <- regsubsets(taste ~ ., data = cheddar[train, ], method = "exhaustive")
summary(model.exh)
## Subset selection object
## Call: regsubsets.formula(taste ~ ., data = cheddar[train, ], method = "exhaustive")
## 3 Variables (and intercept)
## Forced in Forced out
## Acetic FALSE FALSE
## H2S FALSE FALSE
## Lactic FALSE FALSE
## 1 subsets of each size up to 3
## Selection Algorithm: exhaustive
## Acetic H2S Lactic
## 1 ( 1 ) " " " " "*"
## 2 ( 1 ) " " "*" "*"
## 3 ( 1 ) "*" "*" "*"
Definimos…
predict.regsubsets <- function(object, newdata, id, ...) {
form <- as.formula(object$call[[2]])
mat <- model.matrix(form, newdata)
coefi <- coef(object, id = id)
xvar <- names(coefi)
mat[, xvar] %*% coefi
}
Calculamos…
val.errors <- rep(NA, 3)
Y <- cheddar[test, ]$taste
for (i in 1:3) {
Yhat <- predict.regsubsets(model.exh, newdata = cheddar[test, ], id = i)
val.errors[i] <- mean((Y - Yhat)^2)
}
val.errors
## [1] 250.9134 142.4993 177.7624
¿Qué deducimos de aquí?
coef(model.exh, which.min(val.errors))
## (Intercept) H2S Lactic
## -31.121999 3.054151 25.210179
Lactic es sin duda nuestro predictor mas influyente.
Ahora minimizamos el error.
regfit.best <- regsubsets(taste ~ ., cheddar[-obs.out, ])
coef(regfit.best, which.min(val.errors))
## (Intercept) H2S Lactic
## -28.348798 3.315346 22.020523
Ligero cambio
n <- nrow(cheese)
k <- n # numero de grupos, como es elemento a elemento hay n
folds <- sample(x = 1:k, size = nrow(cheese), replace = FALSE)
cv.errors <- matrix(NA, k, 3, dimnames = list(NULL, paste(1:3)))
for (j in 1:k) {
best.fit <- regsubsets(taste ~ ., data = cheese[folds != j, ]) # cogemos datos del train
for (i in 1:3) {
pred <- predict.regsubsets(best.fit, newdata = cheese[folds == j, ], id = i) # datos test
cv.errors[j, i] <- mean((cheese$taste[folds == j] - pred)^2)
}
}
mean.cv.errors <- apply(cv.errors, 2, mean)
Ahora recogemos nuestros términos
mean.cv.errors
## 1 2 3
## 120.27050 68.20906 66.91742
coef(best.fit, which.min(mean.cv.errors))
## (Intercept) Acetic H2S Lactic
## -6.449428 -6.642362 3.884110 29.898710
n <- nrow(cheese)
k <- 4 # numero de grupos
folds <- sample(x = 1:k, size = nrow(cheese), replace = TRUE)
cv.errors <- matrix(NA, k, 3, dimnames = list(NULL, paste(1:3)))
for (j in 1:k) {
best.fit <- regsubsets(taste ~ ., data = cheese[folds != j, ]) # cogemos datos del train
for (i in 1:3) {
pred <- predict.regsubsets(best.fit, newdata = cheese[folds == j, ], id = i) # datos test
cv.errors[j, i] <- mean((cheese$taste[folds == j] - pred)^2)
}
}
mean.cv.errors <- apply(cv.errors, 2, mean)
Comprobamos los términos
mean.cv.errors
## 1 2 3
## 106.09030 66.22008 64.28153
coef(best.fit, which.min(mean.cv.errors))
## (Intercept) Acetic H2S Lactic
## -13.310537 -3.386748 4.131931 20.510849
model.cv <- lm(taste ~ H2S + Lactic, data = cheese)
summary(model.cv)
##
## Call:
## lm(formula = taste ~ H2S + Lactic, data = cheese)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.6439 -4.2257 -0.8544 5.7354 14.7477
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -28.349 8.042 -3.525 0.00191 **
## H2S 3.315 1.114 2.976 0.00697 **
## Lactic 22.021 7.850 2.805 0.01031 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.907 on 22 degrees of freedom
## Multiple R-squared: 0.7395, Adjusted R-squared: 0.7158
## F-statistic: 31.22 on 2 and 22 DF, p-value: 3.751e-07
plot(lm(taste ~ H2S + Lactic, data = cheese), which = 1)
Se observa que practicamente se ajusta por completo a un modelo lineal.
plot(lm(taste ~ H2S + Lactic, data = cheese), which = 2)
Los residuos se comportan como una distribución normal de manera bastante clara
residualPlot(model.cv)
### 5. Errores hay que moverlos aquí
Después de las comparaciones realizadas con respecto a nuestras comprobaciones cruzadas concluimos que el mejor modelo es taste ~ Lactic.
model.final <- model.final1
Es un modelo especialmente simple, sus coeficientes son:
coeff <- summary(model.final)$coeff[,1]
Y algo más de información, aunque ya la hemos estudiado antes la aporta
sumary(model.final)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -39.8244 11.0385 -3.6078 0.001875
## Lactic 42.9502 7.4102 5.7961 1.388e-05
##
## n = 21, p = 2, Residual SE = 9.87699, R-Squared = 0.64
Una consecuencia interesante de tener solo 3 predictores es que podemos interpretar nuestro modelo de regresión como un plano en el espacio.
Primero veamos como se distribuye la variable taste en función de H2S y Lacti
plot_ly(x=H2S, y=Lactic, z=taste, type="scatter3d", color=taste, mode='markers') %>%
layout(scene = list(xaxis = list(title = 'H2S (%)'),
yaxis = list(title = 'Lactic (%)'),
zaxis = list(title = 'Taste (0-100)')))
Y ahora podemos añadir nuestro “plano” de regresión. Marcamos en rojo las observaciones que peor se ajustan al modelo.
# planereg <- scatterplot3d(x=H2S, y=Lactic, z=taste, pch=16, cex.lab=1,
# highlight.3d=TRUE, type="h", xlab='H2S (%)',
# ylab='Lactic (%)', zlab='Taste (0-100)')
# plot.new(planereg$plane3d(model.final, lty.box="solid", col='mediumblue'))
# Me da error